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Abstract 

We consider a system of q diffusing particle species Ai, A2, ■ ■ ■ , Aq that are all equivalent under a 
symmetry operation. Pairs of particles may annihilate according to Ai + Aj —* with reaction rates 
kij that respect the symmetry, and without self-annihilation (ku = 0). In spatial dimensions d > 2 
mean-field theory predicts that the total particle density decays as p{t) ~ provided the system 
remains spatially uniform. We determine the conditions on the matrix k under which there exists a 
critical segregation dimension dseg below which this uniformity condition is violated; the symmetry 
between the species is then locally broken. We argue that in those cases the density decay slows down 
to p{t) ^ f-'^/'^nag £qj- 2 < d < dsog- We show that when dsog exists, its value can be expressed in terms 
of the ratio of the smallest to the largest eigenvalue of k. The existence of a conservation law (as in 
the special two-species annihilation A + B ^ 0), although sufficient for segregation, is shown not to 
be a necessary condition for this phenomenon to occur. We work out specific examples and present 
Monte Carlo simulations compatible with our analytical results. 

PACS 05.40.-a, 82.20.-w 
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1 Introduction 



In nonequilibrium statistical mechanics there exist two annihilation reactions that have 
become prominent model systems. These are the pair annihilation processes A + A ^ 
and A + B ^ 0. Interest in fluctuation effects on the flrst reaction process dates back to 
the seminal work by Smoluchowski p], and in more recent times at least to Bramson and 
Griffeath and for the second process fluctuations were flrst studied by Ovchinnikov and 
Zeldovich The flnal states of both processes are trivial: All particles (save perhaps a 
single one) disappear in the case of A + A — > 0, and only the excess number of particles of 
the initial majority species is left in the case of A + B ^ 0. The nature of the approach 
to the flnal state is, however, nontrivial and has been the subject of many investigations. 

The A + A —>■ problem applies, for example, to domain walls in a one-dimensional 
system like in the zero-temperature Ising model Its asymptotic behavior is the same as 
that of the diffusion-limited coagulation process A + A ^ A . The ensuing asymptotic 
density decay ~ has been experimentally observed (in an intermediate time window) 
in the fusion kinetics of laser-induced excitons in quasi one-dimensional N(CH3)4MnCl3 
(TMMC) polymer chains [B^. 

The study of the two-species annihilation reaction A + B ^ was motivated originally 
by the cosmological problem of matter-antimatter annihilation. Ovchinnikov and Zeldovich 
, and independently a few years later Toussaint and Wilczek jTj , asked whether in such a 
simple annihilation model it would be possible that locally in space only, say, matter would 
be left. The answer rapidly turned out to be positive: The A + B ^ system exhibits the 
phenomenon of species "segregation" , that is, the emergence of ever growing single-species 
domains (either A or B). As a result, after a short initial period, annihilation takes place 
only in "reaction zones" where the domains border, and the decay of the total density 
is slowed down markedly. Indeed, a non-classical decay ~ has been experimentally 
conflrmed in three dimensions in a calcium / fluorophore reaction system [Hj. The condition 
for segregation to occur is that the spatial dimension d be less than the critical segregation 
dimension dseg, which in the case of A + B ^ is equal to ciseg = 4. Early work on 
this two-species problem was done by Kang and Redner clever heuristic reasoning and 
numerical work on the spatial structure of the domains is due to Leyvraz and Redner JU] ; 
and rigorous results were derived by Bramson and Lebowitz jllj . 

The renormalization group approach to the reaction A+A ^ was pursued by Peliti [Sj 
and Lee [12], and to the process A + B Ohy Lee and Cardy [IB], and Oerding [T3j. The 
strong interest in these two reactions has sparked off research on many related reaction- 
diffusion problems by various techniques. Here we will refer only to general introductions 
and overviews, e.g. Refs. [T3]-[?F]. 

In 1986 Ben-Avraham and Redner crediting Kang, introduced a g-species pair 
annihilation problem that interpolates between these two models, and hence puts them in a 
common perspective. They considered a system of q distinct particle species Ai, A2, . . ., Ag, 
all propagating with the same difffusion constant, and reacting according to 

A + Aj^O (1 < 2 < J < g) (1) 

with a single flxed rate. For q = 2 and q = 00 this g-species Mutual Annihilation Model 
(g-MAM) reduces to the well-understood paradigmatic cases A + B and A + A ^ 
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0, respectively. (For the latter case, note that the probability of identical species to 
meet becomes vanishingly small in the limit of g — * cxd; hence any particle encounter will 
lead to a reaction, and the distinction of different species becomes meaningless.) The 
general situation has recently been studied analytically and by Monte Carlo simulations 
in Refs. |2iEniEIl- 

For equal initial densities of all species the asymptotic decay with time t of the total 
particle density p{t) follows a power law ~ t~°'. In Table I we summarize the values of 
the decay exponent a, known or believed to be exact, as a function of q and of the spatial 
dimension d [201 IHIIl IHI] • We also indicate whether or not the density decay is accompanied 
by particle species segregation. The results for c? = 1 and q = 3,4, .. . are based on the 
exact solution of a modified model argued by the authors of Refs. [23 ISIj to be equivalent 
to the one-dimensional g-MAM in the large-time limit. 





q = 2 
A + B ^0 


g = 3,4,... 


q — oo 
A + A^O 


d > 4 


1 


1 


1 


2 < d < 4 


d/4 {*) 


d=2 


1 (log) 


1 (log) 


l<d<2 




d/2 


d= 1 


1/4 (*) 


{q-l)/2q {*) 


1/2 



Table 1. The density exponent a of the asymptotic decay law p{t) ~ t^" for the total partical density in the 
g-species Mutual Annihilation Model, as a function of the number of species q and the spatial dimension d. An 
asterisk (*) indicates that segregation occurs and (log) signifies the appearance of logarithmic corrections at the 
upper critical dimension dc = 2. 

The g-species MAM is characterized by an upper critical dimension dc = 2, below which 
mean-field theory breaks down and renormalization is needed. Physically, in dimensions 
d < 2 the recurrence properties of random walkers lead to depletion zones about each 
surviving particle, and to particle ant icorrelat ions induced by the annihilation kinetics. 
A systematic renormalization group treatment, wherein the spatial dimension d can be 
treated as a continuous variable, has only been carried through for the cases g = oo 
(equivalent to the reaction A + A — > 0) [5J, and for g = 2 (the two-species reaction 
A + B ^ {)) ^3 which explains the empty entries in the center of Table I. 

The derivation of the decay laws listed in Table I relies heavily on the full permutational 
symmetry of the particle species in the g-MAM. It is therefore natural to ask what happens 
to these asymptotic results when the symmetry is lowered. A case of lower symmetry 
occurs, for example, when the particle species are ordered cyclically and annihilation, with 
rate fci, is possible only between particles of two neighboring species along the cycle. The 
reaction constants are then given by 

kij = ki [S(i-j)modq,l + %-i)modg,g-l] (1 < "^j j < g) • (2) 

We will refer to this particular system as the g-species Cyclic Annihilation Model (g- 
CAM). Other examples that display lower than permutational symmetry may easily be 
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constructed. 

In this work we develop, within the mean-field framework {i.e., in dimensions d > 2), 
a general method for finding the exponent of the total density decay law p{t) ~ t~", and 
determining whether or not segregation occurs in a g-species model where the reaction 
rates kij lower the full permutational symmetry. 

In Sec. |21we present our model, which is fully defined by the qxq matrix k of annihilation 
rates. In Sec. El we study how the particle numbers of the individual species fluctuate 
around their average densities at any given time. In Sec. |3 we show, by comparing the 
average decay law with the fluctuations around it, how the reactions come to an end in a 
flnite (i-dimensional volume L"^, and we determine the particle densities in the flnal state. 
In Sec. El we treat the full d- dimensional space as built up from blocks of size L'^ that do 
not communicate for times less than the diffusion time, which is of order L^. This allows 
us to derive our general results concerning the occurrence of segregation and the exponents 
in the density decay law. In Sec. IHl we consider a speciflc example which, by means of a 
suitable parameter A, interpolates between the g-CAM (for A = 0) and the g-MAM (for 
A = 1). In Sec. Owe present Monte Carlo simulations which, although preliminary and 
of limited statistical accuracy, are compatible with our analytical flndings. In Sec. El we 
comment on an alternative approach, namely via rate equations with an additional particle 
diffusion term, before we conclude in Sec. El 

Our considerations allow us to answer another question as well: In the two-species 
A + B ^ system the appearance of segregation is usually explained heuristically through 
arguments invoking the local conservation of the difference between the number of A and 
B particles. This conservation law thus seemed to be at the origin of the segregation 
phenomenon. However, the discovery [21] that the one-dimensional MAM exhibits segre- 
gation for all q < oo, even though it is not subject to any conservation law, cast doubt 
upon this apparent direct link between conservation laws and segregation. This doubt 
was subsequently reenforced by a heuristic argument according to which the MAM should 
exhibit segregation in the entire region of the q-d plane delimited by 1 < g < 3 and 
2 < d < A/{q — 1) jnij. (The only weakness of this argument is that the region concerned 
by it does not contain any points with integer "physical" q and d, except for q = 2, where 
the conservation law applies.) This present work clarifles this issue: We show that the 
existence of a conservation law, although a sufficient condition for segregation in dimen- 
sions d < 4, does not constitute a necessary condition: Segregation may in fact occur in 
the absence of any conservation law. 

2 A general pair annihilation system 

We consider a system of q distinct diffusing particle species Ai, A2, . . . ,Aq subject to the 
pair annihilation reactions 



Here the kij = kji represent the reaction rates per unit of density; they are constrained 
only by the requirements that they be non-negative and that there exist a symmetry 
operation under which all particle species are equivalent. We will set ku = 0, implying 





(3) 
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that no self-annihilation is possible. Without loss of generality we may consider the matrix 
k to be irreducible, i.e., the system of reactants cannot be decomposed into mutually 
noninteracting subsystems. The g-species MAM and the g-species CAM discussed in the 
Introduction are obtained as special cases when all reaction rates in Eq. (jS)) are equal and 
when they are given by Eq. respectively. 

We shall therefore proceed to analyze the mean-field behavior of the system (jS)) for a 
general non-negative traceless symmetric matrix k. Our focus will be on the power law for 
the density decay, p(t) ~ and on the segregation properties of this system. Straight- 
forward scaling analysis tells us that for pair annihilation processes such as considered 
here, very generally mean-field theory is applicable in spatial dimensions d > dc = 2, and, 
likely with logarithmic corrections, at the upper critical dimension dc = 2. We shall oper- 
ate entirely at the mean-field level, and hence our results will concern dimensions d > 2. 
In our discussion we will briefly touch upon implications for lower dimensions d < 2. 

Note that the system of reactions © may be visualized as a graph Q with q vertices 
representing the q species, in which bonds carrying weights kij connect the vertices of pairs 
of species that may react with each other. At certain points in our discussion it will be 
convenient to refer to this graph representation. 

3 Mean-field theory and fiuctuations 

3.1 Averages 

Let the stochastic variable Aj(t) denote the particle number of species Ai present at time 
t in a system of volume L'^, and let riiit) = {Ni(t)) / L'^ be the average density of the for 
i = 1,2, . . . ,q. Here (. . .) is an average with respect to (i) the initial distribution of the iVj 
at time t = and (ii) all realizations of the stochastic time evolution. In order to express 
the mean-field rate equations drii/dt = — kijUiUj in terms of dimensionless variables, 
we set ko = kij (which because of the symmetry between all species is independent of i) 
and no = q~^ Si^i(O). We furthermore define the quantities pi(t) = ni(t)/no, = kij/ko, 
and r = riokot. The mean-field equations then become 

dpi v-^ 

^ = - 2^ i^ij PiPj , (4) 
j 

now with the normalization Ylj '^ij = 1- If initial densities ^^(O) = uq and hence 
Pi{0) = 1 for all i, then the solution of Eq. (jH) reads 

A(r)=p(r) = -^, (5) 
1 + r 

which tends to zero as a power law with decay exponent a = 1. For unequal initial densi- 
ties the asymptotic decay will generically be exponential (readily generalized to stretched 
exponential for c? < 2 as a consequence of the reaction rate renormalization [SlEll2Zj); 
with one or several of the species tending toward a positive limit density pi(cxo); this case 
will not be considered any further below. 
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3.2 Fluctuations 

In principle, the fluctuations around this mean-fleld average are encoded in the master 
equation for the probability distribution P{Ni, N2, . . . , iV^; t), 

E ji[iN. + l){N, + l)Pi...,N, + l,N, + l,...)-N,N,P{...,N,,N,...)]. (6) 

l<i<j<q 

A convenient means to extract them proceeds via van Kampen's Q-expansion j^, as 
was demonstrated for the g-MAM by Ben-Avraham and Redner j2H], and more recently 
for a cyclic three-species trapping reaction by Ben-Naim and Krapivsky [32], and for a 
zero-dimensional population dynamics model by Newman and McKane 

Anticipating that the fluctuations should be of the order of the square root of the total 
particle number Yl'^=i 

N,{t) = noL'^pir) + {n^L'^Y/^ ^,{t) , (7) 

where, due to previous deflnitions, (7i(0)) = 0, and transform the probability distribution 
P on the extensive variables iVj into an equivalent one, to be called F, on the intensive 
variables 7,, according to 

P(iVi, A^a, . . . , N,- 1) = {noLy/^ F{j,, 72, ... , 7^; r) . (8) 

Expanding the master equation (jH)) for P to second order in powers of (noL'^)^^^'^, as in 
Refs. [28 , [33], and and exploiting the rate equations (jl} for pj(r) = p(r) then yields 
a Fokker-Planck equation with time-dependent coefficients, 

drF= ^^^J[p{^^ + ^,){^, + 'yJ) + '^p\d, + d,y]F , (9) 

l<«<j<'7 

where di = d/d'ji. This equation is valid provided the second term on the r.h.s. of Eq. (|7j) 
remains much smaller than the first one. 

3.3 Equations for the second moments 

From the Fokker-Planck equation we may readily deduce equations for the time evolu- 
tion of the averages (7i), the variances (7^^) = Tu, and the covariances {'ji'jj) = Tij. Since 
d(7i)/dr = —pi'ji) — pJ2j given that the (74) are zero initially, we see that 

they vanish for all times. The covariance matrix F satisfies 

— = {Kij + 5ij) - 2p Tij - p^^Kii Tji - p E '^i^ • 

Now let U be the real unitary matrix with elements f/^j that renders k = UkU~^ diagonal, 
and denote its eigenvalues by k^. In the same manner, we define F = UTU~^. Upon 
transforming Eq. (fTUI) to these new variables we obtain 

'^=[l + l{k, + k,)]{p'5,,-2pt^,) . (11) 
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Let us suppose that the particle numbers initially have independent and identical fluctu- 
ations, i.e., Tij{0) = ToSij. Then also Tfj,u{0) = Fq^^^, and it follows from Eq. (|TT|l that 
^ij.u{t) = for yU 7^ z/ at all r > 0. The remaining equation of interest is 

^ = (1 + k,)p' - 2(1 + k,)p T,, . (12) 

Notice that since p depends on time, it is not useful to scale the time r in ()12|1 with the 
factor 1 + K^. The normalization relation Kij = 1 implies that = 1 is an eigenvalue 
of K, with eigenvector (1, 1, ... , 1). By the Perron- Frobenius theorem, and since the are 
necessarily real, we have — 1 < < 1 for all p. In case = —1, the r.h.s. of Eq. ()12|1 
vanishes and the initial fluctuation f^^ does not decay, i.e., an eigenvalue —1 signals the 
presence of a conservation law. 

3.4 Solution of the moment equations 

The solution of Eq. (fT^ . except when = — |, is a sum of two power laws, 

f,,(r) = (Fo-ir,)(l + r)-2(i+'^'') + K,{l + T)-' , (13) 

where we have introduced the coefficient = (l + K^)/(l + 2/t^). Note that the amplitude 
of (only) the first term in Eq. ()13|1 depends on the initial fiuctuation strength Fg. Its 
exponent 2(1 + /t^) varies for the different modes p, which shows that this problem is 
characterized by a spectrum of power laws. The second term in Eq. ()13p is just proportional 
to the density p(r) and hence decays as r~^. Therefore, if < — the first term decays 
more slowly than the second one (and vice versa for > — |). When it so happens that 
= — |, then Eq. (fT!?jl should be replaced with the special logarithmic solution 

T,,{t) = Fo(l + r)-' + i(l + r)-i log(l + r) . (14) 
The variance of the particle number for species i is given by 

i i /I 

Let us denote the algebraically smallest eigenvalue of nhj k^, its degeneracy by Cm, and 
the corresponding value of the coefficient K^j^ by K^. For asymptotically large times we 
then deduce from Eqs. (fT^ - (fT^ the behavior 

ir-Mogr (^kra = -^, r^oo), (16) 

c^q-\To-Kj r-2(i+M 

with the abbreviation K = q^^ ^ i^^. 
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4 Final state in a finite volume 



We are now ready to exploit the results of the preceding sections. We begin by posing the 
question whether there exists a characteristic time r*(L) = noA;ot*(-Z^) at which the particle 
number iVj in the volume L"^ displays a root mean-square deviation from the average that 
is of the order of the average itself, i.e., 



Yet the Fokker-Planck equation was derived under the hypothesis of small fluctuations 
(the second term in Eq. (jZj) was supposed to small compared to the first one), and hence 
ceases to be valid for r ^ t^{L). The issue then is what happens at and beyond this time 
scale. Fluctuations satisfying (fTTjl indicate that the different particles species have widely 
varying numbers. As we argued in the Introduction, in this situation the annihilation 
processes come to an end epxonentially, on the prevailing time scale, through the successive 
extinction of one or several particle species; this final decay is no longer described by the 
Fokker-Planck equation Q. In the finite volume L'^, the final state therefore consists of 
a collection of particles that are not subject to pair annihilations anymore. If all Kij > 
for i ^ j, then at most a single species can survive. But if some of the rates Kij vanish, a 
broader variety of final states is possible (we shall return to this point in SeclH)). In any of 
these cases, whereas initially all species were equivalent, the final state has this symmetry 
broken. In fact, this symmetry breaking may be traced back to the first term in Eq. ()13p. 
which includes the amplitude Fq — K^^. Thus, this symmetry breaking is enhanced by 
initial fluctuations, as represented by Fq, yet it actually persists even for the case Fq = 0, 
i.e. if we initially take all particle numbers Ni{0) to be strictly equal. We shall discuss 
several examples of such broken symmetry states in Sec. IHl 

We turn now to the determination of the characteristic time 'r^,(L). Utilizing in Eq. p7|) 
that {N,{t)) ^ uqL'^t-^ for r ^ oo as well as inserting the asymptotic results of Eqs. p5|) 
and ()16|) . and then solving for t^{L), we find 



where C = [cmQ ^(Fq — -R'm)]^'''^'^™'- The order of magnitude N^{L) of the particle number 
for a surviving species will be N^{L) ~ noL'^p{T^,{L)) ^ hqL'^ / t^,{L) , and hence is given by 



(17) 




(18) 




When Km > — I this surviving number is of order unity, but when = it grows 
logarithmically with the volume; and when < it grows with a positive power of L. 



8 



5 Segregation in an infinite volume 



We now consider an infinite volume in which the q annihilating species propagate with 
(uniform) diffusion constant D, and apply a heuristic analysis: We imagine this infinite 
volume divided into hypercubic subvolumes of size L''. For times less than the characteris- 
tic diffusion time, i.e., for r^rdig(L) ~ L?' / D, the subvolumes can be treated as effectively 
independent, and hence the results of Eqs. flTH|) and (fTIJ|) apply to each of them. We will 
discuss the four distinct cases separately. 

(i) Case — | < < 1- From Eq. (fT^ we infer in this situation that in an isolated vol- 
ume the particle number fluctuations become of the order of the average itself only at times 
when the total particle numbers N^{L) ~ K have decreased to order unity themselves. In 
any case, in dimensions d > 2 the time scale at which this happens, r*(L) ~ L'^, is much 
larger than the diffusion time rdifT(L) ~ L^. Hence in an infinite system the subvolumes 
begin to mix diffusively long before the fluctuations of their particle numbers become of 
the order of the averages. This is tantamount to stating that for /^m > ~| there is no 
particle segregation. 

(ii) Case —l<k^< In this case Eq. fITHj) states that the fluctuations become of 
the order of the average value p(r) at a characteristic time r^=(L) ~ i^d/\2Km\^ g^^^ 
isolated volume the annihilation processes will come to an end on this same time scale. 
Similarly in an infinite system, subject to the condition that 

< rdiff(L) , (20) 

each subvolume of size L'^ will reach a quasi-final state, to be referred to as a "domain", 
and further reactions are possible only at the "reaction zones" separating the domains. 
The quasi-final state in each subvolume is reached independently of its neighbors, and 
the emergence of these different disjoint quasi-final states constitutes the phenomenon of 
segregation into domains. The condition (j20|) implies that segregation occurs in dimensions 
below a critical segregation dimension dgeg 

dseg = 4 l^ml • (21) 

It therefore appears that in the present case the critical segregation dimension for diffusion- 
limited multi-species pair annihilation processes is always constrained by the interval 2 < 
dseg < 4. The dimensionless particle density in a given domain will be of the order p*{L) = 
N ^{L) / {hqL'^) and cannot decrease further until diffusion between neighboring domains 
permits new reactions to take place. Hence we conclude that at any time r the particle 
density in an infinite system equals p*(i^difr(^)), where Ldifr(t) ~ [DtY^"^ = {Dt /nokoy^"^ . 
By combining these relations we find that the density decays asymptotically as 

p(r) ~ (CV)-'^/'^'- {2<d<dseg), (22) 

where C = [cin5'~^(ro — K^)]^'^^'^nQ^'^^^ D / kg is a dimensionless constant. Thus we obtain 
a = d/dseg for 2 < < d^cg- 

(iii) Case Km = — |. In this case dseg = 2, which is also the upper critical dimension dc 
where mean-field theory is only marginally applicable. Nevertheless, upon following the 
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same reasoning as above we find the asymptotic density decay law 

p(r) ~ {k^/2D) log(Dr/A;o) {d = 4eg = 4 = 2) (23) 

in two dimensions. Based on renormalization group arguments and Monte Carlo simula- 
tions, this result was predicted in particular for the 3-MAM |29| IHU} 1^. for which indeed 
Subsec. 16.11 below). It is remarkable that the logarithmic corrections ()23|) 
take the same form as quite generically predicted for pair annihilation reactions in two 
dimensions as consequence of the required reaction rate renormalizations induced by the 
appearance of depletion zones in low dimensions. Although segregation effects are very dif- 
ficult to capture within the framework of renormalized field theory (compare the discussion 
in Ref. |SI]), we therefore hypothesize, but cannot prove, that no additional logarithmic 
factors beyond those exhibited in Eq. appear in the general case. 

(iv) Case = —1. This special case, for which ciseg = 4, occurs when the graph 
Q is bipartite, i.e., when it is comprised of two subsets of vertices such that all bonds 
are between vertices of different subsets. Since these subsets are necessarily equivalent 
under symmetry, this obviously requires that q be even, and we may relabel the sets to 
{1, 2, . . . , q/2} and {g/2 + 1, q/2 + 2, . . . ,q}. The particle numbers then obey the local 
conservation law 

q/2 q 

^Ni{t) - Ni{t) = constant , (24) 

i=l i=q/2+l 

and it is readily verified that k, has an eigenvector (1,...,!,— 1,...,— 1) with eigenvalue 
— 1. For Km = —1 we moreover infer from Eq. p3j) that r^^(r) = Fq, which establishes 
that as a direct consequence of the conservation law the initial fiuctuations in the total 
particle numbers of the two subsets do not relax. Examples where this happens are the 
2- MAM (the two-species pair annihilation reaction A + B ^ 0) and the g-CAM with even 
q (see Subsec. 16. 2|) . 

6 An explicit example 

We proceed to investigate a specific g-species pair annihilation system that depends on a 
parameter A, interpolating for < A < 1 between the g-CAM and the g-MAM that were 
defined in the Introduction. 

We imagine the g species arranged in a cycle. The annihilation rate is set equal to ki 
for any pair of species that are nearest neighbors along the cycle, and to /c2 for any other 
pair. Then the ratio A = /c2/^i is the only intervening parameter in the normalized matrix 
K. Using the abbreviation Afgx = 1 + |(g — 3)A, one has 

K-ij = [A + (1 - A) (5(i_j)modg,l + %-i)modg,g-l)]/2A/'gA (^ 7^ j) (25) 

and = 0. For < A < 1 this model is symmetric only under cyclic translation of the 
species, and the special case A = yields the g-CAM. For A = 1 the model displays full 
permutation symmetry and we recover the g-MAM. 
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It is straightforward to diagonalize k by Fourier transforming in the variable i ~ j, 
and hence the mode labels i/, . . . belong to the set of wavenumbers p G {27in/q} with 
n = 0,l,...,g — 1. The corresponding eigenvalues are 

«:p = g"'5Z/ti,e'P^*-^)= [|A(g5po-l) + {1 - X) cos p] / Af,x ■ (26) 

We shall first discuss the two limiting models, the MAM and the CAM, whose segregation 
properties turn out to be quite different. Next, our study as a function of A will charac- 
terize the crossover between these two cases, thereby shedding light on the nature of the 
segregation mechanism. 

6.1 Mutual Annihilation Model (MAM) 

The parameter value A = 1 defines the Mutual Annihilation Model (g-MAM). In this case 
the graph Q with q vertices is fully connected and every bond carries the same weight. 
This model, introduced in Ref. j^, has been extensively studied in Refs. j^Hl EUl EI] • We 
may reproduce the known results in dimensions c? > 2 as follows. Eq. ()2(j|l shows that for 
A = 1 the matrix k has a single eigenvalue 1 (namely for the mode with p = 0) and a 
(g — l)-fold degenerate eigenvalue Km = —l/{q — 1) (which occurs for all p ^ 0). In the 
two-species model with q = 2, this latter eigenvalue equals —1, whence we encounter the 
special case (iv) of Sec. El There is a conservation law, we find dg^g = 4, and for d < dgeg 
the density decays as p(r) ~ r""^/^ PHUCSIEII- 

For 2 < g < 3, the eigenvalue = —^/iQ — 1) lies in the interval [—1, — |], whence 
we are confronted with case (ii) of Sec. According to Eq. (j2II), there then exists a 
g-dependent segregation dimension 

d,,g{q) = 4/(g - 1) , (27) 

a relation first derived in Ref. .3L. Although this case does not encompass any integer 
values of g, it does suggest that segregation may occur even in the absence of a conservation 
law. Such a scenario will indeed be confirmed below. 

The marginal value g = 3 corresponds to the special case (iii) of Sec. El Indeed, it 
was concluded in Refs. |^ |M] on the basis of renormalization group arguments that in 
two dimensions the total density in the 3- MAM decays as p(r) ~ logr, in agreement 
with Eq. (jSni)- For g > 3 we obtain dseg{q) < 2, and our present theory is not apphcable 
anymore. The analytical and numerical studies reported in Refs. (201 1^ however 
yield species segregation in one dimension, with a g-dependent density decay exponent 
«(g) = (g-l)/2g. 

6.2 Cyclic Annihilation Model (CAM) 

For A = the cyclic graph Q only contains nearest-neighbor bonds, i.e., a particle of 
a given species can annihilate only with particles of one of its two neighboring species. 
This model has been called the Cyclic Annihilation Model or g-CAM. In this case the 
eigenvalues given in Eq. reduce to kp = cos p. For even g, the mode with p = tt yields 
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the smallest eigenvalue = —1. But for odd q the smallest eigenvalue is acquired for 
p = 7i{l ± q~^), namely km = — cos ^. It then follows from Eq. that 

f 4 (g = 2,4,6,...) 

dse, = { . . (28) 
[ 4 cos I (g = 3,5,7,...) 

For g = 2, 3, 4 these results for dseg{q) were previously known, since in these cases the CAM 
actually coincides with a MAM: the 3-CAM coincides with the 3-MAM, and the 2-CAM 
and 4-CAM are both equivalent to the 2- MAM [SHI 1^- The first novel case is therefore 
g = 5, for which Eq. gives 

4eg(5) = 1 + ^5 = 3.236... . (29) 

It follows that the 5-CAM exhibits segregation in spatial dimensions d = 3 and d = 2, and 
Eq. (1221) implies that p(r) ~ r"" with 

^fi(V5- 1)^0.618... (<i^2) 
1 |(V5- 1)= 0.927... (<i = 3) 

(Curiously, the decay exponent in two dimensions is just the golden mean.) In simulations, 
aside from the observation of species segregation, direct measurement of the exponent 
values a provides a straightforward means to verify the present theory (see Sec. 13 below). 

The domain structure in the case of the g-CAM needs to be discussed. Whereas the 
MAM is necessarily characterized by single- species domains, this is no longer true for 
the CAM. Let us first consider g = 5. For the 5-CAM one can certainly conceive the 
possibihty of single-species domains. Any given particle species, however, does not react 
with two other species, and hence is able to coexist with either of those. It appears 
evident, therefore, that a single-species domain is not stable against penetration by either 
of the two species with which its particles cannot annihilate. Thus we must suppose that 
a typical segregated domain will always contain two species, say i and i', out of the five, 
viz. any of the cychcally equivalent combinations {ii'} = {13}, {24}, {35}, {41}, {52}. 
These domains are "stable against penetration" , in the following sense: Any particle of a 
different species intruding into such a domain would annihilate with at least one of the 
two domain species. 

For g > 6 new questions appear that we will not fully address here. As for g = 5, we 
may list the subsets of particles that can coexist in a domain where no interactions take 
place anymore. There are two subsets with three species, {135} and {246}, and three 
subsets with only two species, namely {14}, {25}, and {36}. All these five subsets are 
"impenetrable" to an intruder species. It is therefore an interesting issue which type of 
domains are going to be formed in the segregation process. 

One important conclusion from the g-CAM is that for g = 5, 7, ... we encounter "phys- 
ical" situations (g and d are integers), where, in spite of the absence of a conservation 
law, species segregation into domains emerges. Note that the segregation mechanism does 
not even require initial fluctuations in the species numbers: Even if initially absent, such 
fluctuations are dynamically generated by the annihilation reactions, and they decay more 
slowly than the average particle numbers themselves. 
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6.3 Interpolating model 



We now allow an arbitrary value for A in the interval [0,1], i.e., we consider the more 
general interpolating model. For even q, Eq. tells us that the eigenvalue kp assumes 
its minimum for p = vr; for odd q one may verify that, just as for A = 1, it becomes 
minimal for p as close as possible to vr, that is, for p = 7r(l ±g~^). Upon substituting these 
minimal eigenvalues in Eq. (j2T|) . we obtain 



For A = 0, 1, this expression reduces to the limiting values found above. It is consistent 
with our mean-field assumptions only as long as it leads to a (iseg(?) > 2. That is easily 
seen to imply that at given A, it is meaningful only for q below a maximum value qd^)', 
for q above that value, there is no segregation in any d> 2. 

7 Simulations 

We have performed Monte Carlo simulations in one, two, and three dimensions for the 
g-CAM with q = 2,3,4,5 with the goal to test the theoretical predictions for the decay 
exponent a. For more particle species, our computing resources unfortunately cannot 
provide sufficiently reliable data statistics. 

The simulation algorithm proceeds as described in Ref. jSI]: Starting from an initially 
random distribution of particles on a rf-dimensional {d = 1, 2, 3) hypercubic lattice with 
periodic boundary conditions, the system is evolved as follows: A particle is randomly 
picked. Next, one of its 2d nearest-neighbor sites is selected randomly; if it is occupied by 
a particle of a different species, both particles are removed, otherwise the particle hops to 
the empty site. The Monte Carlo time is scaled with the total number of particles present 
(asynchronous time update). 

Beginning with the data in three dimensions, obtained from simulations on a 100 x 
100 X 100 cubic lattice by averaging over 20 runs with random initial conditions, we plot 
the effective decay exponent 



in Fig. ^ In d = 3, naturally our statistics is worst, and the data cease to be reliable at 
around t ~ 600. The results for even and odd q are clearly distinct; and as anticipated, 
the processes for q = 2 and g = 4 are equivalent. Whereas Eqs. (j^H|) and (f^^ predict 
a = 3/4 asymptotically for q even, the data show that the effective exponent a{t) rather 
reaches a plateau at ~ 0.88. As discussed in Ref. [HI], one should however expect the 
asymptotic decay law to be somewhat masked by the mean- field behavior ~ t~^, which 
could explain the combined effective decay exponent closer to 1. Indeed, we found the 
same deviation from the asymptotic value 0.75 in our simulations for the 2-MAM (see 
Fig. 8 in Ref. jHI])- The data for odd q are similarly plagued by crossover effects. For 
g = 3, for which the MAM and CAM are equivalent, we see a slow convergence towards 




q even 



q odd 



(31) 
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Figure 1: Effective decay exponent a{t) for the three-dimensional g-CAM with q = 2, 3, 4, and 5 
(random initial conditions, equal initial particle numbers for each species). The data were obtained on a 
100 X 100 X 100 cubic lattice, averaged over 20 runs. 



a = 1. The data for g = 5 might indicate that a{t) reaches a plateau at a smaller value, 
in line with the prediction a ~ 0.93 of Eq. (fHUjl . but the deteriorating statistics preclude 
a definite conclusion. 

Our two-dimensional results, from 20 runs on a 1000 x 1000 square lattice, are depicted 
in Fig. 121 The data become unreliable at about t ~ 1000. Once more, the g-CAM for 
q = 2 and g = 4 are seen to be equivalent, with a{t) settling towards the asymptotic value 
1/2, albeit masked again by the competing mean-field power law, just as for the 2- MAM 
(see Fig. 6 in Ref. [SI])- Yet now the runs for q = 3 and q = 5 yield manifestly different 
power laws. In the three-species CAM, the effective exponent is still changing in the time 
window accessible to our simulations, running towards a = 1, perhaps with logarithmic 
corrections as predicted by Eq. (j^ . For g = 5, however, we find a plateau value ~ 0.71, 
perhaps with a slowly decreasing tendency. This may be interpreted as a combination of 
the predicted asymptotic decay law (jHUI) with the mean-field result ~ t^^. 

The simulation data in d = 1, from averaging over again 20 runs on 10^ lattice sites, 
are reliable at least up to t ~ 10"^, and again confirm the equivalence of the g-CAM 
with even g. For both g = 2 and g = 4 we obtain a slow but definite approach towards 
the predicted a = 1/4, as for the 2-MAM (see Fig. 2 in Ref. jSI]). For the 3-CAM, 
equivalent to the 3-MAM, we find the expected slow convergence towards a = 1/3 from 
above fl^ EOl 1^- Remarkably, the data for g = 5, for which no analytical prediction 
is available in one dimension, appear to reach a{t) ~ 0.33 as well, but faster and from 
below. This is remarkably close to what Eq. ()22|) would predict if we applied it (without 
justification) in dimension d = 1, namely a = l/(igeg(5) = 0.309..., with the segregation 
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Figure 2: Effective decay exponent a{t) for the two-dimensional g-CAM with g = 2, 3, 4, and 5 (random 
initial conditions, equal initial particle numbers for each species). The data were obtained on a 1000 x 1000 
square lattice, averaged over 20 runs. 
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Figure 3: Effective decay exponent a{t) for the one-dimensional g-CAM with g = 2, 3, 4, and 5 (random 
initial conditions, equal initial particle numbers for each species). The data were obtained on a 10^ lattice, 
averaged over 20 runs. 
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dimension given in Eq. (jSHI)- 

In summary, our simulation results for the g-CAM, within the accuracy and statistical 
errors of our data, are consistent with the analytical predictions of the previous sections, 
if we allow for crossover effects induced by the presence of the competing mean-field decay 
law. The equivalence of the 4-CAM and 2-CAM, which are both identical to the 2-MAM 
(A + i? — > 0), is clearly demonstrated. Since the 3-CAM and 3-MAM are identical, only our 
data for g = 5 yield novel results. Unambiguous identification of the density exponents a 
would ideally require better statistics, that is, simulations on considerably larger systems. 

8 Differential equation approach 

Before concluding we wish to make a methodological remark. One may attempt to address 
the issue of segregation by considering space-dependent particle densities pj(x, r) that 
satisfy the react ion- diffusion equations (with diffusivity D) 

— — — = DApi - K,j piPj . (33) 

j 

For spatially uniform (well-mixed) systems these reduce to the ordinary differential equa- 
tions (jH). Let us now apply to the system (jH^ what is essentially a linear stability analysis. 

To this end, we write pi{x,T) = p(t)[1 + ei(a;,r)] and linearize with respect to the 
perturbations ej. With Eqs. (jlj) this yields 

^^^^ = DAe.ix, r) - p(r) ^ K,,e, , (34) 

j 

where the neglected nonlinear terms are of order pe^. In terms of the linearly transformed 
variables 

e,{p, r) = J2 I ^'x e'P-' e,(x, r) (35) 

with U as defined in Subsec. 13. 3| the time evolution equation becomes diagonal and 
reads 



[Dp' + p{r)k^]e^{p,T) . (36) 



For given initial e^(p, 0) this is solved by 

i,{p,T) = e,{p,0){l + T)-'^e-''P'^ . (37) 

Irrespective of the value of this solution exists and tends to zero as r ^ oo. Nevertheless, 
for < —Dp^ it will increase initially. In the event that |e^| grows with time such as to 
be no longer negligible with respect to unity at some instant of time, this signals that the 
neglected nonlinear terms in the differential equation begin to play a role; in particular, 
they will prevent from decreasing beyond —1 and hence from turning negative. We 
will argue, however, that in this case the physical justification for the differential equations 
()33|) breaks down, and that in fact segregation sets in. 
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This goal requires going beyond the mathematics and invoking the interpretation of 
the density as an average over discrete particles. For particles initially randomly and 
uniformly distributed in a volume L'^, we typically have that e^(p, 0) ~ L~'^/'^. If this initial 
fluctuation is negative, and if in the course of time e^{p,T) approaches —1, this should 
be interpreted properly as there being an appreciable probability for species extinction. 
A necessary condition for the r.h.s. of Eq. (jHTj) to become of order unity reads explicitly 
L~'^/^(l + T)^'^^' > 1, or, since must be negative, r > = L"'/'^"'''!. But on these 
time scales Dp'^r must still be much less than unity. Taking the smallest allowed valued 
for p in the volume L'^, that is, ~ 27r/L, the additional condition Dp^r^ <^ 1 yields 
j^d/\2Kf,\-2 ^ which will happen when d < 4|k^|. 

Hence the differential equation approach, combined with appropriate considerations in 
which the particle discreteness intervenes, leads to the identical expression for the critical 
segregation dimension as found in Eq. ()2H) of Sec. El The differential equation approach 
was followed, essentially, in Ref . . To our opinion, however, the Fokker-Planck method 
used in Refs. j2Hl ES] and in the present work is preferable, since it is based directly on 
the more fundamental description of an interacting many-body problem in terms of a 
master equation. More specifically, the stochasticity taken into account in the differential 
equation approach is due only to the random fluctuations present in the initial state. As a 
consequence, this approach wrongly suggests that the initially dominant species (or, more 
precisely, their initially dominant mode fi) is also the surviving one. 

It follows from the Fokker-Planck equation that there is actually an interplay between 
initial and dynamically generated fluctuations, brought out by the solution (fT^ of the 
second moment equations. This solution contains a dynamically generated contribution to 
the fluctuations [viz. the terms oc K^), and a contribution due to the initial fluctuations 
(the Fq term). In systems without segregation {i.e., when > — |) the dynamically 
generated fluctuations become larger than those due to the initial conditions, which are 
eventually forgotten. In systems with segregation (for < — |) the initial and dynamical 
contributions are of the same order and the initial conditions co-determine the final state. 
In this weaker sense, the segregation phenomenon that appears in the cases with < 
is still linked to the persistence of initial fluctuations, even in the absence of conservation 
laws. 

9 Conclusion 

We have studied a wide class of g-species reaction-diffusion systems with pair annihilation 
processes between distinct species. This class includes the two well-studied paradigmatic 
cases A + A ^ {in the limit q oo) and A + B ^ {q = 2). Within mean- field 
theory, i.e. for spatial dimensions d > 2, we have determined for each member of the class 
(i) whether or not segregation occurs, and (ii) the value of the decay exponent a in the 
asymptotic power law for the total particle density. Our findings represent a considerable 
extension of previous work on segregation in diffusion- limited annihilation reactions. Our 
preliminary simulation data are compatible with the analytical results. 

Our method builds on ideas that were applied earlier in the context of g-species models 
by Ben-Avraham and Redner |2HI, and more recently also by Ben-Naim and Krapivsky 
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[nni and Newman and McKane This approach is obviously not hmited to the special 
cases studied here, but we have not aimed at being exhaustive. It is straightforward, for 
example, to drop the restriction of no self-annihilation {ka = 0), and redo the analysis. 
The method may also be employed to analyze more involved situations, such as models 
with species belonging to two distinct equivalence classes. 

Various open problems remain, in particular concerning the nature of the domain struc- 
ture in several important cases. Furthermore, the theory presented here is unable to ad- 
dress the issue of what happens in low dimensions, d < 2, where particle anticorrelations 
become manifest and render the mean-field treatment invalid. We expect that more elabo- 
rate simulations will be carried out in the future to test our predictions as well as to study 
questions beyond reach of the present theory. 
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